#### "Inclusionary regimes, party institutionalization and redistribution under authoritarianism" ####
# authors: "Lars Pelke"
# date: 2020-06-16
# written under "R version 3.6.1 (2019-07-05)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

#libraries
library(tidyverse)
library(ggpubr)
library(texreg)
library(dotwhisker)
library(broom.mixed)
library(countrycode)

library(lme4)
library(sjPlot)
library(sjmisc)
library(lfe)

# My ggplot2 theme
theme_lp <- theme(
  legend.position = "bottom",
  panel.background = element_rect(fill = NA),
  # panel.border = element_rect(fill = NA, color = "grey75"),
  axis.ticks = element_line(color = "grey95", size = 0.3),
  panel.grid.major = element_line(color = "grey95", size = 0.3),
  panel.grid.minor = element_line(color = "grey95", size = 0.3),
  legend.key = element_blank())

# set working directory
# please use the working directory, where you stored the zip-file. 

#### Load Data ####
vdem <- readRDS("data/vdem_merged.rds") 
vdem_spaw <- readRDS("data/vdem_merged_spaw.rds") 

vdem <- vdem %>%
  drop_na(v2x_regime)

vdem_spaw <- vdem_spaw %>%
  drop_na(v2x_regime)

#### Generate Sample ####

vdem <- vdem %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_spaw <- vdem_spaw %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_rel <- vdem %>%
  drop_na(rel_red)

vdem_spaw <- vdem_spaw %>%
  drop_na(universalism_all)

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE) # Distinct observations
vdem_spaw <- distinct(vdem_spaw, country_id, year, .keep_all= TRUE) # Distinct observations

vdem_rel <- vdem_rel %>% # delete those countries with less than 5 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_spaw <- vdem_spaw %>% # delete those countries with less than 5 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_rel <- vdem_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_spaw <- vdem_spaw %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_rel2 <- vdem_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_spaw2 <- vdem_spaw %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_rel <- vdem_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, rel_red, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_rel2 <- vdem_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, rel_red, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

vdem_spaw <- vdem_spaw %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, universalism_all, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_spaw2 <- vdem_spaw2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, universalism_all, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

#### Computing group-mean and de-meaned variables ####

vdem_rel <- sjmisc::de_mean(vdem_rel, pol_incl, v2xps_party, v2x_regime,
                            v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                            grp = country_id)

vdem_rel2 <- sjmisc::de_mean(vdem_rel2, pol_incl, v2xps_party, v2x_regime,
                             e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

vdem_spaw <- sjmisc::de_mean(vdem_spaw,pol_incl, v2xps_party, v2x_regime,
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

vdem_spaw2 <- sjmisc::de_mean(vdem_spaw2, pol_incl, v2xps_party, v2x_regime,
                              e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                              v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                              grp = country_id)

vdem_rel$auto_regime_type <- factor(vdem_rel$auto_regime_type, 
                          levels = c(0, 1, 2),
                          labels = c("Closed Autocracy", 
                                     "Hegemonic Multiparty Autocracy", 
                                     "Competitive Multiparty Autocracy"))

vdem_spaw$auto_regime_type <- factor(vdem_spaw$auto_regime_type, 
                                    levels = c(0, 1, 2),
                                    labels = c("Closed Autocracy", 
                                               "Hegemonic Multiparty Autocracy", 
                                               "Competitive Multiparty Autocracy"))


#### Table Summary Authoritarian Regime Types ####

#### Table A 6 and Table A7 ####

## Redistribution Dataset ##
vdem_rel_regime_type <- vdem_rel %>%
  group_by(auto_regime_type) %>%
  count()

vdem_rel_regime_type <- vdem_rel %>%
  group_by(auto_regime_type) %>%
  summarize(n = n_distinct(country_name))


vdem_rel_regime_type2 <- vdem_rel2 %>%
  group_by(auto_regime_type) %>%
  count()

vdem_rel_regime_type2 <- vdem_rel2 %>%
  group_by(auto_regime_type) %>%
  summarize(n = n_distinct(country_name))

## Universalism Index ##

vdem_spaw_regime_type <- vdem_spaw %>%
  group_by(auto_regime_type) %>%
  count()

vdem_spaw_regime_type <- vdem_spaw %>%
  group_by(auto_regime_type) %>%
  summarize(n = n_distinct(country_name))


vdem_spaw_regime_type2 <- vdem_spaw2 %>%
  group_by(auto_regime_type) %>%
  count()

vdem_spaw_regime_type2 <- vdem_spaw2 %>%
  group_by(auto_regime_type) %>%
  summarize(n = n_distinct(country_name))

#### Figure A3 Density rel redistribution ####

ggplot(data = vdem_rel) +
  geom_density(aes(x = rel_red), fill = "grey60") +
  facet_wrap(~auto_regime_type) +
  theme_classic() +
  labs(y = "Density",
       x = "Relative Redistribution",
       title = "") 
ggsave("Output/A3.pdf", height = 10, width = 20, units= c("cm"))

#### Figure A4 Density universalism index ####

ggplot(data = vdem_spaw2) +
  geom_density(aes(x = universalism_all), fill = "grey60") +
  facet_wrap(~auto_regime_type) +
  theme_classic() +
  labs(y = "Density",
       x = "Universalism Index",
       title = "") 
ggsave("Output/A4.pdf", height = 10, width = 20, units= c("cm"))

#### Figure A5 Density inclusion ####

ggplot(data = vdem_rel) +
  geom_density(aes(x = pol_incl), fill = "grey60") +
  facet_wrap(~auto_regime_type) +
  theme_classic() +
  labs(y = "Density",
       x = "Political Inclusion by regime type",
       title = "") 
ggsave("Output/A5.pdf", height = 10, width = 20, units= c("cm"))

#### Figure A6 Density party insitutionalization ####

ggplot(data = vdem_rel) +
  geom_density(aes(x = v2xps_party), fill = "grey60") +
  facet_wrap(~auto_regime_type) +
  theme_classic() +
  labs(y = "Density",
       x = "Party Institutionalization",
       title = "") 
ggsave("output/A6.pdf", height = 10, width = 20, units= c("cm"))

#### Figure A8 Inclusion, Party Institionalization and Redistribution ####

plot1 <- ggplot(data = vdem_rel, aes(x = pol_incl, y = v2xps_party)) +
  geom_jitter(alpha = 0.3) +
  geom_smooth(method = lm, colour = "black") +
  theme_bw() + 
  labs(y = "Party Institutionalization",
       x = "Political power by socio-economic and socital group",
       title = "") 

plot2 <- ggplot(data = vdem_rel, aes(x = pol_incl, y = rel_red)) +
  geom_jitter(alpha = 0.3) +
  geom_smooth(method = lm, colour = "black") +
  theme_bw() + 
  labs(y = "Relative redistribution",
       x = "Political Inclusion",
       title = "") 
ggarrange(plot1, plot2, ncol = 1, nrow = 2)

plot3 <- ggplot(data = vdem_rel, aes(x = v2xps_party, y = rel_red)) +
  geom_jitter(alpha = 0.3) +
  geom_smooth(method = lm, colour = "black") +
  theme_bw() + 
  labs(y = "Relative redistribution",
       x = "Party Institutionalization",
       title = "")

ggarrange(plot1, plot2, plot3, ncol = 1, nrow = 3)

ggsave("Output/A8.pdf", height = 20, width = 15, units= c("cm"))


#####################################################################################
#####################################################################################

#### Sample description "Redistribution Dataset", Countries and Years ####

sample1 <-vdem_rel2 %>% 
  group_by(country_name) %>% 
  summarize(start = min(year), 
            end = max(year), 
            no_year = n())

sample1$start <- as.character(sample1$start)
sample1$end <- as.character(sample1$end)

cols1 <- c("country_name", "start", "end", "no_year")

library(stargazer)
stargazer(as.data.frame(sample1[, cols1]), type = "latex", summary = FALSE)

#### Sample description "Universalism Dataset", Countries and Years ####

sample2 <-vdem_spaw2 %>% 
  group_by(country_name) %>% 
  summarize(start = min(year), 
            end = max(year), 
            no_year = n())

sample2$start <- as.character(sample1$start)
sample2$end <- as.character(sample1$end)

cols1 <- c("country_name", "start", "end", "no_year")

library(stargazer)
stargazer(as.data.frame(sample2[, cols1]), type = "latex", summary = FALSE)

#### Figure A2 ####

ggplot(data = vdem_rel, aes(x = pol_incl, y = v2xps_party)) +
  geom_jitter(alpha = 0.3) +
  geom_smooth(method = lm, colour = "black") +
  facet_wrap(~auto_regime_type) +
  theme_classic() +
  labs(y = "Party Institutionalization",
       x = "Political inclusion",
       title = "") 
ggsave("Output/A2.pdf", height = 10, width = 20, units= c("cm"))

#### Figure A7 ####

vdem_rel$country_year <- str_c(vdem_rel$country_name, vdem_rel$year, sep = " ")
  
library(ggrepel)

ggplot(data = vdem_rel, aes(x = pol_incl, y = rel_red)) +
  geom_jitter(alpha = 0.3) +
  geom_label_repel(aes(label=ifelse((pol_incl < -1 | pol_incl >2 ) & (rel_red <= 0 | rel_red >= 30)  ,as.character(country_year),'')),hjust=0,vjust=0) +
  geom_smooth(method = lm, colour = "red") +
  facet_wrap(~auto_regime_type) +
  theme_classic() +
  labs(y = "Relative redistribution",
       x = "Political Inclusion",
       title = "") 
ggsave("Output/A7.pdf", height = 15, width = 20, units= c("cm"), dpi = 1200)


#### Figure A1: Maps with Redistribution and inclusion / Party institutionalization ####

library(ggplot2)
library(dplyr)
require(maps)
require(viridis)
theme_set(
  theme_void()
)


#### Load Data ####
vdem <- read_csv2("calculations/data/vdem_merged.csv") 

## Load World Map ##

world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill="white", colour = "black")

# Change data for rel_red, party institutionalization and rel red 

vdem_rel <- vdem 

vdem_rel <- vdem_rel %>%
  mutate(rel_red = ifelse(v2x_regime<=1, rel_red, NA),
         pol_incl = ifelse(v2x_regime<=1, pol_incl, NA), 
         v2xps_party = ifelse(v2x_regime<=1, v2xps_party, NA))

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE) # Distinct observations

vdem_rel <- vdem_rel %>% # delete those countries with less than 5 country-year observations to deal with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 3) %>%
  ungroup(country_id)


#Reduce data to mean_data_frame
vdem_map_rel <- vdem_rel %>%
  filter(year>=2000) %>%
  group_by(country_name) %>%
  summarize(pol_incl_m = mean(pol_incl, na.rm = TRUE), 
            rel_red_m = mean(rel_red, na.rm = TRUE),
            v2xps_party_m = mean (v2xps_party, na.rm = TRUE)) %>%
  rename(region = country_name) %>%
  mutate(region = ifelse(region == "United States of America", "USA", region)) 
  

#Merge data with maps

vdem.map <- left_join(vdem_map_rel, world_map, by = "region")

# Mapping 
title <- theme(
  plot.title = element_text(hjust = 0.5),
  plot.subtitle = element_text(hjust = 0.5), 
  legend.position = "right"
)

#Relative redistribution
map_red <- ggplot(vdem.map, aes(long, lat, group = group))+
  geom_polygon(aes(fill = rel_red_m), color = "black")+
  scale_fill_viridis_c(option = "inferno", na.value="white") +
  labs(fill = "",
       title = "Relative Redistribution in Autocracies",
       subtitle = "averaged between 2000 and 2018") + 
  title 


# Inclusion
map_incl <- ggplot(vdem.map, aes(long, lat, group = group))+
  geom_polygon(aes(fill = pol_incl_m), color = "black")+
  scale_fill_viridis_c(option = "inferno", na.value="white") +
  labs(fill = "",
       title = "Political Inclusion in Autocracies",
       subtitle = "averaged between 2000 and 2018") + 
  title 

# Party Institutionalization
map_party <- ggplot(vdem.map, aes(long, lat, group = group))+
  geom_polygon(aes(fill = v2xps_party_m), color = "black")+
  scale_fill_viridis_c(option = "inferno", na.value="white") +
  labs(fill = "",
       title = "Party Institutionalization in Autocracies",
       subtitle = "averaged between 2000 and 2018") + 
  title 

ggarrange(map_red, map_incl, map_party, ncol = 1, nrow = 3)
ggsave("output/A1.pdf", height = 30, width = 20, units= c("cm"))


  
